%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This Matlab program is for 'The explicit formula for the Hodrick-Prescott filter in finite sample' by Adriana Cornea-Madeira
%This probram compares the results from Corollaries 6 and 7 with the HP filter based on the direct inversion 
%of the matrix in (2)
clear all; clc;

n = 11; %sample size 
l = 2;
dates = (1:1:n)';

c = randn(n,1); 
tau = zeros(n,1); eta = randn(n,1);
tau(1) = eta(1);
for i = 2:n
    tau(i) = tau(i-1) + eta(i);
end
Y = tau + c; %data or use load('data')

%smoothing parameter
alpha = 5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HP filter weights computed using our formulae from Corollaries 6 and 7, 
% for sample size even and odd
m = floor(n/2);
i = linspace(1,m,m)';
j = i;
T1 = sqrt(2/(n+1))*sin(i*j'*pi/(n+1));%matrix of eigenvectors
x = repmat(i',m,1);
J = ones(m,m);
J = J.*(-1).^(x-1); 
  
M2 = sqrt(2/(n+1))*sin(i*(2*j)'*pi/(n+1)); 
G2 = (2/(n+1))*((2*sin((2*j)*pi/(n+1))-sin(2*(2*j)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j)*pi/(n+1))).^2))*((2*sin((2*j)*pi/(n+1))-sin(2*(2*j)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j)*pi/(n+1))).^2))';
k2 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*sin((2*j)*pi/(n+1))-sin(2*(2*j)*pi/(n+1))).^2./(1+alpha*(2-2*cos((2*j)*pi/(n+1))).^2)));
E = k2*M2*G2*M2';

% eigenvalues
b = cos(j*pi/(n+1));
gam1 = 1./(1+4*alpha*(1-b).*(1-b));
p = m:-1:1;
Pm = eye(m);
Pm = Pm(p,:);
gam2 = 1./(1+4*alpha*(1+Pm*b).*(1+Pm*b)); %if n odd, the (m+1)'s eigenval is left out but we do not need it here.

W1 = T1*diag(gam1)*T1 + (T1.*J)'*Pm*diag(gam2)*Pm*T1.*J;
W2 = T1*diag(gam1)*(T1.*J)'*Pm+(-1)^l*(T1.*J)'*Pm*diag(gam2)*Pm*(T1.*J)'*Pm.*J;
   
if m == n/2
    % if n is even
    V1 = W1;
    V2 = W2;
    j = i-1;
    k1 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1))).^2./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2)));
    M1 = sqrt(2/(n+1))*sin(i*(2*j+1)'*pi/(n+1)); 
    G1 = (2/(n+1))*((2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2))*((2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2))';
    D = k1*M1*G1*M1';
    H = V1+D+E;
    % the inverse
    Finv_new = [H V2+(D-E)*Pm;Pm*(V2*Pm+D'-E') Pm*H*Pm]; 
    % trend estimate
    tau_est = [H*Y(1:m)+(V2+(D-E)*Pm)*Y(m+1:end);Pm*(V2*Pm+D'-E')*Y(1:m)+Pm*H*Pm*Y(m+1:end)];
else
    % if n is odd
    xm1 = sqrt(2/(n+1))*sin(i*(m+1)*pi/(n+1));
    xm2 = sqrt(2/(n+1))*sin((i+m+1)*(m+1)*pi/(n+1));
    lam11 = 1/(1+alpha*(2-2*cos((m+1)*pi/(n+1)))^2);
    V1 = W1 + lam11*xm1*xm1';
    V2 = W2 + lam11*xm1*xm2';
    i = linspace(1,m,m)';
    j = [i-1;m];
    k1 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1))).^2./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2)));
    M1 = sqrt(2/(n+1))*sin(i*(2*j+1)'*pi/(n+1)); 
    G1 = (2/(n+1))*((2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2))*((2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2))';
    D = k1*M1*G1*M1';
    s = i; t = j; j = m+1; 
    d = zeros(m,1); e = zeros(m,1);
    for i = 1:m
        d(i,1) = k1*sum(sum(((2/(n+1))^2)*((2*sin((2*t+1)*pi/(n+1))-sin(2*(2*t+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*t+1)*pi/(n+1))).^2).*sin(i*(2*t+1)*pi/(n+1)))*((2*sin((2*t+1)*pi/(n+1))-sin(2*(2*t+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*t+1)*pi/(n+1))).^2).*sin(j*(2*t+1)*pi/(n+1)))'));
        e(i,1) = k2*sum(sum(((2/(n+1))^2)*((2*sin((2*s)*pi/(n+1))-sin(2*(2*s)*pi/(n+1)))./(1+alpha*(2-2*cos((2*s)*pi/(n+1))).^2).*sin(i*(2*s)*pi/(n+1)))*((2*sin((2*s)*pi/(n+1))-sin(2*(2*s)*pi/(n+1)))./(1+alpha*(2-2*cos((2*s)*pi/(n+1))).^2).*sin(j*(2*s)*pi/(n+1)))'));
    end
    xm1m1 = sqrt(2/(n+1))*sin((m+1)*(m+1)*pi/(n+1));
    v = T1*diag(gam1)*xm1 + xm1*lam11*xm1m1+(T1.*J)'*Pm*diag(gam2)*xm2;
    v_m1 = xm1'*diag(gam1)*xm1 + lam11*xm1m1^2 + xm2'*diag(gam2)*xm2;
    i = m+1; 
    d_m1 = k1*sum(sum(((2/(n+1))^2)*((2*sin((2*t+1)*pi/(n+1))-sin(2*(2*t+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*t+1)*pi/(n+1))).^2).*sin(i*(2*t+1)*pi/(n+1)))*((2*sin((2*t+1)*pi/(n+1))-sin(2*(2*t+1)*pi/(n+1)))./(1+alpha*(2-2*cos((2*t+1)*pi/(n+1))).^2).*sin(j*(2*t+1)*pi/(n+1)))'));
    e_m1 = k2*sum(sum(((2/(n+1))^2)*((2*sin((2*s)*pi/(n+1))-sin(2*(2*s)*pi/(n+1)))./(1+alpha*(2-2*cos((2*s)*pi/(n+1))).^2).*sin(i*(2*s)*pi/(n+1)))*((2*sin((2*s)*pi/(n+1))-sin(2*(2*s)*pi/(n+1)))./(1+alpha*(2-2*cos((2*s)*pi/(n+1))).^2).*sin(j*(2*s)*pi/(n+1)))')); 
    H = V1+D+E;
    % the inverse
    Finv_new = [H (v+d+e) V2+(D-E)*Pm;(v+d+e)' (v_m1+d_m1+e_m1) (v+d-e)'*Pm;Pm*(V2*Pm+D'-E') Pm*(v+d-e) Pm*H*Pm]; 
    % trend estimate
    tau_est = [H*Y(1:m)+(v+d+e)*Y(m+1,1)+(V2+(D-E)*Pm)*Y(m+2:end);(v+d+e)'*Y(1:m,1)+(v_m1+d_m1+e_m1)*Y(m+1,1)+(v+d-e)'*Pm*Y(m+2:end,1);Pm*(V2*Pm+D'-E')*Y(1:m)+Pm*(v+d-e)*Y(m+1,1)+Pm*H*Pm*Y(m+2:end)]; 
end
c_est = Y - tau_est; % cyclical component

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HP filter weights computed using the HP filter based on the direct
% inversion of the matrix in (2)
p = n:-1:1;
P = eye(n);
P = P(p,:);

Q = diag(2*ones(n,1))+diag(-1*ones(n-1,1),1)+diag(-1*ones(n-1,1),-1);% or use tril&tril triu(ones(4,4),-1)
g = [-2;1;zeros(n-2,1)];
A = eye(n)+alpha*Q*Q;
F = A-alpha*g*g'-alpha*P*g*g'*P;
Finv_old = inv(F); % Finv_old should be the same as Finv
trend = Finv_old*Y; % trend should be the same as tau_est
c_old = Y - trend; % c_old should be the same as c_est
